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Abstract 

The Fluctuation Theorem gives an analytical expression for the probability of observing second 
law violating dynamical fluctuations, in nonequilibrium systems. At equilibrium statistical mechanical 
fluctuations are known to be ensemble dependent. In this paper we generalise the Transient and 
Steady State Fluctuation Theorems to various nonequilibrium dynamical ensembles. The Transient 
and Steady State Fluctuation Theorem for an isokinetic ensemble of isokinetic trajectories is tested 
using nonequilibrium molecular dynamics simulations of shear flow. 
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1. INTRODUCTION 

The fluctuation theoremi-3 (FT) gives a general formula for the logarithm of the probability 
ratio that the time averaged dissipative flux takes on a value, Jj , to minus the value, namely - J; , in a 
nonequilibrium steady state. This formula is an analytic expression that gives the probability, for a 
finite system and for a finite time, that the dissipative flux flows in the reverse direction to that 
required by the Second Law of Thermodynamics. For steady state trajectories, the Steady State FT 
(SSFT) is only true in the long time limit. Evans and Searles^ have shown that if transient trajectories 
are considered rather than steady state trajectories, a Transient FT (TFT) that is true at all times can be 
derived. Further, when the nonequilibrium steady state is unique, one would expect the asymptotic 
convergence of the Transient to the Steady State Fluctuation Theorem since averages over transient 
segments should then approach those taken over nonequilibrium steady state segments. However, 
there has been some recent discussion4,5 of this point and not all parties agree on asymptotic 
convergence. 

The Transient FT considers the thermostatted response of an ensemble of systems to an applied 
dissipative field. The system is thermostatted so that it may reach a steady state after a Maxwell time, 
Xm- In the original derivations, it was supposed that the initial ensemble was the microcanonical 
ensemble and the dynamics was isoenergetic. This initial ensemble was chosen because the 
probability of observing trajectories originating in a specified phase volume is simply proportional to 
the measure of that volume. However, it is straightforward to apply the same procedure to an 
arbitrary initial ensemble, and an analytical form of the FT that is valid at all times can be obtained. In 
this note we demonstrate how the FT can be extended to an arbitrary system and as an example, 
derive a transient FT for the isokinetic nonequilibrium response of an initial isokinetic ensemble. 

2. TRANSIENT FLUCTUATION THEOREM (TFT) 

Consider an N-particle system with coordinates and peculiar momenta, 
{qi,q2,..qN,Pi,"PN} -(Q^P) -T. The internal energy of the system is 

N 

Hq = V pf / 2m + $(q) = K + 4> where $(q) is the interparticle potential energy, which is a function 

i=i 
of the coordinates of all the particles, q and K is the total peculiar kinetic energy. In the presence of 



an external field Fg, the thermostatted equations of motion are taken to be, 

q,=p,/m + Q(r)F, 

p. = F,(q) + D,(r)F,-a(r)p, 

where Fj(q) = -3$(q) / dq^ and a is the thermostat multiplier which in this case is applied to the 
peculiar momenta, and Ci and Di represent the coupling of the system to the field. The dissipative 
flux is given by J where H^'' = -J(r)VFg. 

The probability that a trajectory segment will be observed within an infinitesimal phase space 
volume of size 5V about T at time t, Pr(5V(r(t),t)) is given by, 

Pr(5v(r(t),t)) = f(r(t),t)5v(r(t),t) . (2) 

where f(r(t),t) is the normalised phase space distribution function at the point r(t) at time t. The 
Lagrangian form of the non-equilibrium Kawasaki distribution function2 is given by: 



f(r(t),t) = exp[-£ A(r(s))ds]f(r(0),0) = exp[-£' A(r(s))ds 



f(r(0),o) (3) 



where A(r) = 3f • /dT is the phase-space compression factor. Now consider the set of initial 
phases in the volume element of size 5V(r(0),0) about r(0) . At time t, these phases will occupy a 
volume 5V(r(t),t) . Since the number of ensemble members within a comoving phase volume is 
conserved, (3) can be used to show. 



5V(r(t),t) = expf£'A(r(s))ds 



5V(r(0),0). (4) 



which is simply the phase space volume contraction along the trajectory, from r(0) to r(t) . 

We will refer the to the trajectory starting at r(0) and ending at r(t) as T(^^iy As discussed 
previously2, a time-reversed trajectory segment that is initiated at time zeE, r*(0;t)=MF)(r(0;t)) where 
M(T) is the time-reversal mapping, can be constructed by applying a time-reversal mapping at the 



midpoint of this trajectory and propagating forward and backward in time from this point for a period 
of t/2 in each direction. See [2] for further details. The pdht r*(0) is related to the point r(t)by a 
time-reversal mapping. Since the Jacobian of the time-reversal mapping is unity, the phase volume 
5V(r(t), t) is equal to the phase volume 5V(r * (0), 0) . The ratio of the probability of observing the 
two volume elements at time zero is: 

Pr(5v(r(0), 0)) f (r(0), o)5v(r(0), o) 



Pr(5v(r (0),0)) f(r (0),0)5v(r (0),o) 

(5) 



f(r(0),o) 
f(r(t),o) 



exp 



-£A(r(s))ds 



where we have used the symmetry of the time-reversal mapping and equation (4) to obtain the final 
equality. This TFT is completely general and applies to any ensemble or type of dynamics. If the 
initial phase space distribution function is known (regardless of whether it is an equilibrium 
distribution), we can then obtain an analytical expression for the probability ratio. Note that the phase 
space distribution function in the numerator and denominator both refer to that at time zero, therefore 
it is readily applied to a system which is at equilibrium at time zero but moves away from equilibrium 
(when the distribution function may be intractable). 

For arbitrary initial ensembles and arbitrary dynamics (constant energy, temperature and/or 
pressure etc) it is convenient to define a general dissipation function B(r) , so that 



f ds B(r(s)) = in( ^(^(Q^'QM - rA(r(s))ds 
Jo frrcf^ n^ Jo 



f(r(0),0) 

f(r(t),o) 

(6) 



B,t 



we can obtain a formula for the probability ratio of observing a particular value of B( and its 
negative. This is achieved by summing over all appropriate regions of phase and it is straightforward 
to show that a TFT of the form. 
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l„^li^^ = A. (7) 

Pr(-B, = A) 



for the property B is obtained. 2,7 



The equation (7) has previously been considered for a system that is initially in a microcanonical 

ensemble and which undergoes isoenergetic dynamics. 2,6-8 in this case B(r) = -A(r) = -VFg(3J(r) 

Pr(A = A) 

and simulations have verified the resulting TFT, that is In = = -At = -VF BJ,t where 

^ Pr(-A, = A) '^ ' 

Jj = A / (VFgP) .1-2,5,7,8 The purpose of this paper is to derive equations equivalent to (7) for other 
ensembles and dynamics. 



As an example, we consider a system initially in the the isokinetic ensemble and undergoing 
isokinetic dynamics. The isokinetic distribution function is. 



f(r(0),o) = f,(r(0),o) = exp(-PH.,(r(0)))5(K(r(0))-K,) ^^^ 

jdrexp(-|3H„(r))5(K(r)-Ko) 



Substituting into equation (5) gives 

Pr(5v(r(0),o)) _ f^(r(0),0)5v(r(0),0) 



Pr(5V(r(0),0)) fjr(0),0)5V(r(0),0) 



exp(-pH.,(r(0))) r_p ^ ^^ 



exp(-|3Ho(r(t))) 



(9) 



= exp((3£6(r(s))ds)expf-£A(r(s))ds 



where we have used the symmetry of the mapping, Ho(r*(0)) = Ho(r(t)), V(r*(0),0) = V(r(t),t) 

and K(r*(0),0) = K(r(t),t) to obtain the second equality and that 

Ho(r(t)) = Ho(r(0)) + f Ho(r(s))ds = Ho(r(0)) + f'6(r(s))ds to obtain the final equality. We see 

JO *'0 



that: 

B(r) = |36(r)-A(r) 

= -|3J(r)VF, 
and from (7) we therefore have, 

In P^lJt=A) ^_AtF[3v. (10) 

Pr(J, =-A) 

The TFT given by equation (10) is true at all times for the isokinetic ensemble when all initial phases 
are sampled from an equilibrium isokinetic ensemble. 

If the system reaches a unique steady state, then at long times the value of J will fluctuate about 
its steady state value. A set of nonequilibrium time averaged currents, {Jj ^^}, can be generated by 
evolving time along a single phase space trajectory which starts at some initial phase which is 
consistent with the macroscopic conditions (N,V,T or E, Fe, etc). One waits for a time which is 
much longer than the Maxwell time, x,^, which characterises the transient response, before one 
begins to analyse time contiguous "steady state" trajectory segments and computes for example the 
statistics of the set {Jtssl- 

If the steady state is unique, the statistics of the set {J^ g^} is independent of the original initial 
phases. Since the steady state is unique we can also gather the statistical information on {J^ g^} by 
studying the response of an ensemble of initial phases, provided we wait several Maxwell times 
before gathering "steady state" data. After the transient response time, the instantaneous system 
properties will be characteristic of the steady state. The set {i^j} averaged over transient trajectory 
segments, starting from an equilibrium ensemble of phases will approach the (unique) steady state 
set {J, ss) with {J,^} = {J, 55} + 0(x„ / t). 

The transient flux, Jj t ~ ^t ss + ^('^m / t), can be expanded using a Taylor series analysis and 
because at sufficiently long times the deviation of J, jg from its mean value decreases as 0( t '^ ), 
equation (10) can be written, 



In P^^J'T-A) ^^^ Pr(3,ss-A) +iJo(^ /t)t) = -AtFjV. (11) 

Pr(J.,,=-A) Pr(J.,ss=-A) ^ ^ ^ ^^ 



Therefore since Im is finite, provided t is sufficiently long and the steady state is unique we expect, 



lim-ln^S^^^-^^ = -AF^pV . (12) 

<--t Pr(J,ss = -A) 



Although there is some contention regarding this point4,5, the numerical tests presented below 
support the possibility of this convergence. 

In Table I we give Transient FTs for various ergodically consistent ensembles6, that is for 

systems where the zero field dynamics preserves the initial ensemble: df(T,0)/dt Ip ^o~ *-* • ^^ ^^^^ 

give the exact form of the steady state FT derived from the transient FT by assuming the steady state 
is unique. 

The steady state FT, proposed and tested by Evans, Cohen and Morrissi and proven by 
Gallavotti and Cohen referred to particular conditions (that the dynamics is dissipative, isoenergetic,9 
reversible and chaotic) and may be expressed by the formula^ 

ltoll„2£Si^^ = _A. (13) 

t^oc t Pr(At = -A) 

In a Gaussian isoenergetic system the phase space contraction rate is instantaneously proportional to 
both the entropy production rate per unit volume and the dissipative flux, and therefore the 
fluctuations in these three properties will be directly related. Alternative forms of the FT for the 
isoenergetic system in terms of the fluctuations of the dissipative flux or the entropy production rate 
per unit volume are therefore trivially obtained. In the original work of Gallavotti and Cohen, for 
example, the FT was expressed in terms of the entropy production rate per unit volume. 3 For a 
system undergoing isokinetic dynamics the dissipative flux is no longer instantaneously proportional 
to the phase space contraction rate (although the time averaged values are proportional in the limit 



t^oo since A((r) = VFg(3Jt(r) + 0(1 / t)), and it is therefore of interest to consider (13) for the 
isokinetic dynamics. 

3. NUMERICAL RESULTS FOR TRANSIENT AND STEADY STATE SYSTEMS 

We test equations (10) and (13) for systems consisting of N=32 WCA particles in two Cartesian 
dimensions undergoing isokinetic shear flow in two Cartesian dimensions using the SLLOD 
equations of motion for planar Couette flow. 6 The equations of motion are: 



(14) 



qi=Pi/mi+i7yi 
Pi = Fi-iYPyi-api 

where y is the applied strain rate, and a is the Gaussian isokinetic thermostat multiplier: 

N 
E^'i-Pi-'i^xiPyi 

a = ^^^ . (15) 

LPi Pi 

i=l 

In this case the dissipative flux, J is equal to Pxy, the xy element of the pressure tensor and 
A(r) = -2Na(r) + 0(1). All results below are presented in reduced units. The simulations are 
carried out for systems under two conditions. In the first case a temperature of T = 1 .0, a particle 
density of n = 0.8 and an applied strain rate of Fe = Y= 0.5 are employed. This applied strain rate is 
sufficiently high that we are in the nonlinear regime and the shear thinning is approximately 12%. In 
the second case a temperature of T = 1.0, a particle density of n = 0.4 and an applied strain rate of Fg 
= Y = 0.01 are employed and the system is in the linear regime. 

For both sets of conditions, transient trajectories and steady state trajectory segments are 
considered. The transient trajectories of length t are initiated from an isokinetic equilibrium ensemble 
and the time-averaged value of the thermostat multiplier and the dissipative flux is calculated for each 
trajectory. The steady state trajectory segments are initiated at equally spaced time origins along a 
single steady state trajectory. In this case the thermostat multiplier and the dissipative flux are 
calculated for each trajectory segment of length t. 
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Results for transient trajectories are shown in Figures 1 and 2. In Figure 1, transient trajectory 
segments are of length t = 0.6 are considered. The value of Wj(t) = -ln(Pr(Jj)/Pr(-Jj) / (ty|3V) is 
plotted as a function of J, to test equation (10), and W„(t) = ln(Pr(a()/Pr(-aj)) / (2Nt) is plotted as 
a function of OC; to test equation (13). If the equations are valid, a straight line of unit slope is 
obtained. As expected, the results indicate that equation (10) is valid whereas equation (13) is clearly 
not valid for an averaging time t = 0.6. In Figure 2 we show how the slopes of the lines formed 
from a plot of Wj(t) versus J( and W„(t) versus OC; vary with averaging time, t. We denote the 
slopes of these lines Sj(t) and S„(t) respectively. This figure demonstrates that the Transient FT 
derived in this paper for an isokinetic systems and given by (10) is valid at all times, but that equation 
(13) is not valid at the times considered. The behaviour of S„(t) is consistent with a 1 / Vt 
convergence to 1 for the system where Fe = 0.5: this is the same rate at which the standard deviation 
of the distribution of the value of a^ goes to zero, however it can also be fitted to other functional 
forms. The data for the system where Fe =0.01 do not appear to converge to a slope of 1 although it 
has clearly not reached its limiting behaviour at the times considered. 

In Figure 3, we test the FT for steady state simulations where neither equation (10) nor (13) is 
true instantaneously, but we expect (10) to be true for times greater than a several Maxwell times 
(Xm=0.084 for T=1.0, n=0.8, Y=0.5). The steady state trajectory segment averages are obtained from 
a single trajectory, sampled at different points along the trajectory and consider the same state points 
as examined for the transient simulations. Clearly, convergence to the limiting behaviour predicted by 
equation (10) is observed. In contrast, the limiting behaviour predicted by equation (13) has not been 
realised even at the longest times considered, and the data indicate that «/ convergence does occur, it is 
very slow compared to its convergence in an isoenergetic system (note that the n = 0.8 state point is 
quite similar to that studied by Evans, Cohen and Morrissi where convergence to within the error 
bars of the data (3%) was observed when the averaging time was 0.5 and Y= 0.5). 

The values of S„(t) determined in the transient experiments are also shown in Figure 3 to 
demonstrate the convergence of the distribution of the averages over transient trajectories to those 
over steady state trajectory segments with time. Although (13) is not valid at the times considered. 
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the values of S„(t) obtained from the transient experiments agree with those obtained from the steady 
state trajectory segments at sufficiently long times. Furthermore the Sj(t) for the transient and steady 
state segments agree and equation (10) is verified in both cases. These observations support the 
possibility that the steady state relationships are correctly predicted from the limiting behaviour of the 
TFT expressions. It implies that not only does the form of the TFT change when isokinetic dynamics 
replaces isoenergetic dynamics, but that a corresponding new SSFT is obtained which is given by 
equation (12). 

At equilibrium, equation (13) is clearly false for isokinetic dynamics for any finite t. In this case 
the probability of observing positive and negative fluctuations in A^ must be equal and therefore 

S„(t) = for all finite times. Equation (13) incorrectly suggests that the logarithm of the probability 
ratio is proportional to the fluctuation, and that S„(t) = 1 . The falsity of equation (13) at equilibrium 
is consistent with results in figure 3 where the calculated value of S„(t) is less than unity for the times 
considered and the departure is greater for the system which is subject to a smaller field. It seems to 
indicate that either the convergence time becomes infinite as the field goes to zero or that S„(t) 

converges to a value less than unity. The behaviour of the fluctuations in Jj is correctly predicted by 
equation (10) which becomes trivial. The observation that equation (13) is not valid for an isokinetic 
system at equilibrium is not in conflict with Gallavotti and Cohen4 because the equilibrium system is 
both non dissipative and non isoenergetic9 as is required. 

We note that it is difficult to test the FT in nonequilibrium systems for long averaging times, 
large fields or large systems since in these limits, the variance of the distributions approach zero. 
Therefore, the fluctuations become so narrow that by the time the formula has converged to its 
asymptotic behaviour, it is not possible to observe the Second Law violating trajectory segments 
either in a computer simulation or experimentally. This highlights the utility of a Transient rather than 
the Steady State FT. Tests of the asymptotic Steady State FT can of course only be carried out for 
long averaging times. The Transient FT on the other hand, can be tested for arbitrarily short averaging 
times where the probability of spontaneous Second Law violations is much greater, approaching 0.5 
as t^O. 
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The integrated form of the fluctuation formula (IFT) also provides a means of reducing the 
statistical error when testing the FT for long periods of time. 8 The IFT corresponding to equation 
(11) is: 



p(Jt<0) 
P(J.>0) 



exp(J.tFjV))_ ^^ = (exp(J.tFjV)):'^^ (16) 



where p( Jj > 0) and p( Jj < 0) refer to the probabilities of observing trajectory segments with positive 
or negative values of Jj, respectively and the notation (...)j ^^ and (...)j ^^ refer to ensemble averages 
over trajectories with positive or negative values of Jj, respectively. We define 
Yj(t) = (ln(p(Jj > 0) / p(J, < 0)) - 1) / (ln((exp(J,tF^|3V))_ ^^)-l) and in Figure 4 it is verified that 
Yj(t) = 1 at all times, in agreement with (16). (Note that subtraction of 1 in the numerator and 
denominator is carried out to circumvent large statistical errors when the probabilities are almost 
equal. This property therefore differs from Y(t) used to measure convergence of the IFT in reference 
[5]). The time that can be considered in the IFT is still limited, particularly at high fields. 

By defining Y„(t) ^ (ln(p(a. < 0) / p(a. > 0)) - 1) / (ln((exp(-2Na,t))_^„) - 1) an IFT 
corresponding to (13) is 

lim Y„(t) = lim((ln(p(a. < 0) / p(a. > 0)) - 1) / (ln(exp(-2Na.t))^ ^^ - 1)) = 1 (17) 

The behaviour of Y^,(t) for the nonequilibrium systems examined above is shown in Figure 4. 
Convergence to the limiting behaviour indicated by (17) has not occurred at the longest times 
considered. An IFT, valid at all times and for arbitrary field strengths is: 

Pp^ = (exp(J>,pV))^^^^ (18) 

p(at<0) ^ "^.>" 

This relation predicts that Y„(t) = (ln(p(a( < 0) / p(a( > 0)) - 1) / (ln((exp(J,tF^(3V))_ ^g) - 1) is unity 

at all times. Note that in a numerical test of equation (18) the large contribution of rarely observed 
events (when 1^ is highly positive) to the ensemble average is problematic. lo This is observed in the 
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data presented in figure 4. 



4. CONCLUSION 

In this paper we have shown how the Transient FT, which was originally derived for the 
isoenergetic nonequilibrium response of an equilibrium microcanonical ensemble may be generalised 
to cases where the response is observed for different ensembles with different types of thermostats. 
The form of the FT may change with ensemble, as shown in Table 1 . This should come as no 
surprise since even at equilibrium fluctuation formulae are generally ensemble dependent. 

We have presented numerical results for the response of the isokinetic ensemble to isokinetic 
shear flow for a system far from equilibrium in the nonlinear regime and for a system in the linear 
regime. We find that the Transient FT derived here and given by equation (10) for this 
ensemble/dynamics combination is satisfied for all averaging times. From this Transient FT we 
obtain a corresponding Steady State FT given by equation (12) for this system. It is observed that the 
form of the Steady State FT derived from the Transient FT also varies with ensemble. In the limit of 
long averaging times the computational results again support the validity of the proposed Steady State 
FT for this combination of ensemble and dynamics. 
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Figure Captions 

Figure 1. The nomialised logarithmic probability ratio W(X) for the fluctuations in a transient 
response numerical experiment with T = 1.0, n = 0.8, N=32, y= 0.5 and a transient trajectory 
segment of length t = 0.6. The circles are a test of the transient FT given by equation (11) and the 
crosses show that the limiting expression given by equation (13) is not applicable, at least at this value 
of t. The straight line is the expected result from equation (11) and (13). 

Figure 2. The squares show the slope of straight line fitted through a plot of Wj(t) versus Jj, as a 
function of t for the transient response of a 2 dimensional system of 32 particles at a) T = 1.0, n = 0.8 
and Y= 0.5 and b) T = 1.0, n = 0.4 and y= 0.01. Equation (11) predicts a slope of 1 at all times, 
which is indeed observed. The circles show the slope of a straight line fitted through a plot of W„(t) 
versus of a,, as a function of t for the same system. Equation (13) says the slope will be 1 in the 
limit of long times. 0(1/N) effects have been accounted for. 

Figure 3. The squares show the slope of straight line fitted through a plot of Wj(t) versus Jj, as a 
function of t for a steady state simulation of a 2 dimensional system of 32 particles at a) T = 1.0, n = 
0.8 and y = 0.5 and b) T = 1.0, n = 0.4 and y = 0.01. Equation (11) predicts a slope of 1, which is 
consistent with the results. The circles show the slope of a straight line fitted through a plot of W„(t) 
versus of a,, as a function of t for the same system. The crosses show how the slope of a straight 
line fitted through a plot of W„(t) versus a, varies with t using data for a transient experiment (as 
also shown in figure 2). Convergence of the transient and steady state results is observed. 0(1/N) 
effects have been accounted for. 



Figure 4. Tests of the integrated fluctuation formulae given by equations (16)-(18) for a 2 
dimensional system of 32 particles at a) T = 1.0, n = 0.8 and y = 0.5 and b) T = 1.0, n = 0.4 and y = 
0.01. The squares show Yj(t) which is predicted to be unity at all times by equation (16). If 
equation (17) is correct, then Y„(t), shown by the circles, should approach unity at long times. The 
crosses show Y„(t), which will equal unity at all times if (18) is valid. Large numerical errors 
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in Y„(t) are observed for the system at large field, as expected, lo 



Table I. Transient fluctuation formula in various ergodically consistent ensembles.'*'' 



Isokinetic dynamics 
Isobaric-isothermal 
Isoenergetic 

Isoenergetic boundary driven flow 

Nose-Hoover (canonical) dynamics 
Wall ergostatted field driven flow'= 

Wall thermostatted field driven flow'^ 



Steady state isoenergetic dynamics:'' 

1 ft„+t 
J, = - J(s)ds where tg » T^ 
t ^1 



in P^^J^) = 
Pr(-J,) 


-JttFePV 


in P^(Jl) = 
Pr(-Jt) 


-JttFePV 


in '^%^ - 


= -JPttFeV 


, Pr(At) 
or In =!^ = 


-A,t 



Pr(-At) 



j^_Pr(A^_ 
Pr(-At) 



or In — ^^^ = -A,t 
Pr(-At) 

Pr(J.) 
Pr(-J,) 

Pr(Jp.) 



^""P^) = -'^''^^^"^ ~ \{^^A^^'^^ - P^y- I Pw.,)])j^ 



In 



Pr(-Jp.) 
-jp,tF^V-ln 



exp 



F V( ['" J(s)B(s)ds + [""''' J(s)B(s)ds) 

Jo -^tn+t 



Umilnim = -Jp.F^V 
— t Pr(-Jp,) ^' = 



-JFeV = 



d< 
dt 

lyad 

-JFeV=^^ 
dt 

_jFeV = ^ 
dt 



JFeV = 


dH;^'' 
dt 


JFeV = 


dt 


JFeV = 


dt 


JFeV = 


dH;;'' 



dt 



^ It is assumed that the limit of large system has been taken so that 0(1/N) effects can be neglected. In the Nose-Hoover dynamics, it is also assumed that the mass associated with the 

heat bath is an extensive variable. 

^ Hq is the equilibrium internal energy and Iq is the equilibrium enthalpy. 

*= In the wall ergostatted/thermostatted systems, it is assumed that the energy/temperature of the full system (wall and fluid) is fixed. 

'' Similar steady state formula can be obtained for other ensembles. Xy^^ is the Maxwell time that characterises the time required for relaxation of the nonequilibrium system into a steady 

state. 
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